package test;

import org.apache.commons.math.MathException;
import org.apache.commons.math.distribution.NormalDistributionImpl;

public class Bs2cirvolas
{
	private static double bscp=0;
	private static double deviation=1e-7;

	public Bs2cirvolas(double bscp){
		this.bscp=bscp;
	}

	public Bs2cirvolas(double bscp, double deviation){
		this.bscp=bscp;
		this.deviation=deviation;
	}
	
	
	public static void main(String[] args)
	{
		double sigmabs=0.5;
//		 bscp <- pnorm(1/2*sigmabs) - pnorm(-1/2*sigmabs)
		try
		{
			Bs2cirvolas bs2cirvolas = new Bs2cirvolas(sigmabs);
			bscp =pnorm(sigmabs*.5)-pnorm(-sigmabs*.5);
			/**/double d1 = derivateFirst(1);
			double d2 = derivateSecond(1);
//			System.out.println("derivateFirst="+d1);
//			System.out.println("derivateSecond="+d2);
//			
//			System.out.println("f(x,y)= "+ function(3.0,3.0));
//			System.out.println("dz/dx = "+ derivateXFirst(2, 3));
//			System.out.println("dz/dz = "+ derivateYFirst(2, 3));

//			System.out.println("f(x,y)= "+ function(2,3));
//			System.out.println("f(x,y)= "+ function(2,4));
//			System.out.println("f(x,y)= "+ function(2,5));
//			System.out.println("dz/dy = "+ derivateYFirst(2, 3));
			System.out.println("findZero = " +bs2cirvolas.findZero(3,3) );
			
			
			//System.out.println("zero = "+ zeroLine(0.5,2,2));			
			//System.out.println("minimo: "+bs2cirvolas.findMinimum(9));
		}
		catch (MathException e)
		{
			e.printStackTrace();
		}
	}

	private static double derivateFirst(double d) throws MathException
	{
		double prev=Double.MAX_VALUE;
		double res=0;
		for(int i=1;i<10;i++){
			double h = 1/Math.pow(10, i);
			res=(function(d+h)-function(d-h))/(2*h); 
			double rel = Math.abs((prev - res) / (res == 0 ? 1 : res));
			//System.out.println(h+") der="+res+" deviation="+rel+" f(x)="+function(d)+" f(x+h)="+function(d+h));
			System.out.println("\t\t"+i);
			if (rel<deviation)
				break;	
			prev=res;
		}		
		return res;
	}

	private static double derivateXFirst(double x, double y) throws MathException
	{
		double prev=Double.MAX_VALUE;
		double res=0;
		for(int i=1;i<10;i++){
			double h = 1/Math.pow(10, i);
			res=(function(x+h,y)-function(x-h,y))/(2*h); 
			double rel = Math.abs((prev - res) / (res == 0 ? 1 : res));
			System.out.println(h+") der="+res+" deviation="+rel+" f(x)="+function(x-h,y)+" f(x+h)="+function(x+h,y));
//			System.out.println("\t\t"+i);
			if (rel<deviation)
				break;	
			prev=res;
		}		
		return res;
	}

	private static double derivateYFirst(double x, double y) throws MathException
	{
		double prev=Double.MAX_VALUE;
		double res=0;
		for(int i=1;i<10;i++){
			double h = 1/Math.pow(10, i);
			res=(function(x,y+h)-function(x,y-h))/(2*h); 
			double rel = Math.abs((prev - res) / (res == 0 ? 1 : res));
			System.out.println(h+")  f(x-h)="+function(x,y-h)+" f(x+h)="+function(x,y+h)+" diff="+(function(x,y+h)-function(x,y-h)));
//			System.out.println("\t\t"+i);
			if (rel<deviation)
				break;	
			prev=res;
		}		
		return res;
	}
	



	private static double derivateYSecond(double x, double y) throws MathException
	{
		double prev=Double.MAX_VALUE;
		double res=0;
		for(int i=1;i<10;i++){
			double h = 1/Math.pow(10, i);
			res=(function(x,y+h)-2*function(x,y)+function(x,y-h))/(h*h); 
			double rel = Math.abs((prev - res) / (res == 0 ? 1 : res));
			System.out.println(h+") der="+res+" deviation="+rel+" f(x)="+function(x,y-h)+" f(x+h)="+function(x,y+h));
//			System.out.println("\t\t"+i);
			if (rel<deviation)
				break;	
			prev=res;
		}		
		return res;
	}
	


	private static double derivateXSecond(double x, double y) throws MathException
	{
		double prev=Double.MAX_VALUE;
		double res=0;
		for(int i=1;i<10;i++){
			double h = 1/Math.pow(10, i);
			res=(function(x+h,y)-2*function(x,y)+function(x-h,y))/(h*h); 
			double rel = Math.abs((prev - res) / (res == 0 ? 1 : res));
			System.out.println(h+") der="+res+" deviation="+rel+" f(x)="+function(x-h,y)+" f(x+h)="+function(x+h,y));
//			System.out.println("\t\t"+i);
			if (rel<deviation)
				break;	
			prev=res;
		}		
		return res;
	}

	private static double derivateSecond(double d) throws MathException
	{
		double prev=Double.MAX_VALUE;
		double res=0;
		for(int i=1;i<10;i++){
			double h = 1/Math.pow(10, i);
			res=(function(d+h)-2*function(d)+function(d-h))/(h*h); 
			double rel = Math.abs((prev - res) / (res == 0 ? 1 : res));
			//System.out.println(h+") der="+res+" deviation="+rel+" f(x)="+function(d)+" f(x+h)="+function(d+h));
			if (rel<deviation)
				break;	
			prev=res;
		}		
		return res;
	}
	
	
	

	private static double forw(double x, double y) throws MathException{
		double res=0;
//			forw <- function(asw) {
//			 a <- max(asw[1],0);
//			 sw <- abs(asw[2]);
//			  wmin <- sqrt(a) - 1
//			  sw/sqrt(2*pi)*exp(-wmin^2/(2*sw^2))*(wmin + 2) +
//			  (sw^2 + 1-a)*pnorm(-wmin/sw)
//			 }
		double a = Math.max(x, 0);
		double sw = Math.abs(y);
		double wmin = Math.sqrt(a)-1;
		res=sw/Math.sqrt(Math.PI)*Math.exp(-Math.pow(wmin, 2)/(2*sw*sw))*(wmin+2)+(sw*sw+1-a)*pnorm(-wmin/sw);	
		return res;
	}

	private static double pnorm(double d) throws MathException
	{
		NormalDistributionImpl n = new NormalDistributionImpl(0.0, 1);
		return n.cumulativeProbability(d);
	}

	private static double atmcp(double x, double y) throws MathException{
//		 atmcp <- function(asw) {
//		 a <- asw[1]; 
//		 sw <- asw[2];
//		  wmin <- sqrt(1 + a) - 1
//		  sw/sqrt(2*pi)*exp(-wmin^2/(2*sw^2))*(wmin + 2) +
//		  (sw^2-a)*pnorm(-wmin/sw)
//		 }  
		double res=0;
		double a=x;
		double sw=y;
		double wmin=Math.sqrt(1+a)-1;
		res=sw/Math.sqrt(Math.PI)*Math.exp(-Math.pow(wmin, 2)/(2*sw*sw))*(wmin+2)+(sw*sw+1-a)*pnorm(-wmin/sw);
		return res;
	}

	private static double function(double x, double y) throws MathException{
		double res=0;
//		(forw(x)-1)^2 + (atmcp(x)-bscp)^2
		//res=Math.pow(forw(x,y),2)+Math.pow((atmcp(x,y)-bscp), 2);		
		
		
		
		res = x*x+y*y -32;
		res =  2*x*x*x+8*x*x+2*x-12;
		return res;
	}

	private static double function(double x){
		return 2*x*x*x+8*x*x+2*x-12;

		//return 2*x*x+2*x-4;
		//return Math.log(x)+x;
		//return Math.sqrt(x);
	}
	
	private static double zeroLine(double m, double x, double y){
		return (m*x-y)/m;
	}

	private String findZero(double x, double y) throws MathException{
		double rel=0;
		double prev=Double.MAX_VALUE;
		double dx=derivateXFirst(x, y);
		double dy=derivateYFirst(x, y);
		double z=function(x,y);
		double x0=zeroLine(dx, x, z);		
		double y0=zeroLine(dy, y, z);
		double z0=function(x0,y0);
		rel = Math.abs((prev - z0) / (z0 == 0 ? 1 : z0));
		int i=0;
		while(rel>deviation ){
			i++;
			prev=z0;
			dx=x0-z0/derivateXFirst(x, y);
			dy=y0-z0/derivateYFirst(x, y);
			z=function(x,y);
			x0=zeroLine(dx, x, z);		
			y0=zeroLine(dy, y, z);
			z0=function(x0,y0);
			System.out.println(i+") dx="+dx + " dy="+dy+" rel="+rel);
			rel = Math.abs((prev - z0) / (z0 == 0 ? 1 : z0));
			if (i>100)
				break;
		}
		return x0+","+y0;
	}
	
	private double findZero(double x) throws MathException{
		double rel=0;
		double prev=Double.MAX_VALUE;
		double d=derivateFirst(x);
		double fx=function(x);
		
		
		double x0=zeroLine(d, x, fx);
		double y0=function(x0);
		
		rel = Math.abs((prev - y0) / (y0 == 0 ? 1 : y0));
		int i=0;
		while(rel>deviation ){
			i++;
			prev=y0;
			d=derivateFirst(x0);
			fx=function(x0);
			x0=zeroLine(d, x0, fx);
			y0=function(x0);
			rel = Math.abs((prev - y0) / (y0 == 0 ? 1 : y0));
			System.out.println(i+") x="+x0 + " f(x)="+y0+" rel="+rel);
			if (i>100)
				break;
		};
		return x0;		
	}
	
	private double findMinimum(double x) throws MathException{
		double rel=0;
		double prev=Double.MAX_VALUE;
		double d=derivateSecond(x);
		double fx=derivateFirst(x);
		double x0=zeroLine(d, x, fx);
		double y0=derivateFirst(x0);
		rel = Math.abs((prev - y0) / (y0 == 0 ? 1 : y0));
		int i=0;
		while(rel>deviation ){
			i++;
			prev=y0;
			d=derivateSecond(x0);
			fx=derivateFirst(x0);
			x0=zeroLine(d, x0, fx);
			y0=derivateFirst(x0);
			rel = Math.abs((prev - y0) / (y0 == 0 ? 1 : y0));
			System.out.println(i+") x="+x0 + " f(x)="+y0+" rel="+rel);
			if (i>100)
				break;
		};
		return x0;
		
	}
	
}
